##########################################################################################
# Title: Figures — Section 2 (Figure 4 & Figure A1)   
# Section: 2 (Figures)
# Purpose: Faceted bar charts for relative density (to 5 km) and deviations vs. global avg
# IO:
#   IN : Processed_Data/Intermediate/Section_2/Processed_BD_Analysis/ALL_CITIES_100K_09_2020.rds
#   OUT: Graphs/Figure_4.png, Graphs/Figure_A1.png
##########################################################################################

rm(list = ls()); options(stringsAsFactors = FALSE)

suppressPackageStartupMessages({
  library(data.table)
  library(ggplot2)
  library(grid)   # unit()
})

# ---- Project root (align with run_all.R) -------------------------------------
if (!requireNamespace("here", quietly = TRUE)) install.packages("here", quiet = TRUE)
ROOT      <- normalizePath(here::here(), winslash = "/")
IN_DIR    <- file.path(ROOT, "Processed_Data", "Intermediate", "Section_2", "Processed_BD_Analysis")
GRAPH_DIR <- file.path(ROOT, "Graphs")
dir.create(GRAPH_DIR, recursive = TRUE, showWarnings = FALSE)

# ---- Read --------------------------------------------------------------------
infile <- file.path(IN_DIR, "ALL_CITIES_100K_09_2020.rds")
stopifnot(file.exists(infile))
DT <- readRDS(infile) |> as.data.table()

# ---- Prep --------------------------------------------------------------------
DT <- DT[city_name != "Mumbai"]
DT[, bounds_f := factor(bounds, levels = c(5,10,15,20,25,30))]

# (Optional) shorten SA facet labels so they fit
DT[city_name == "Tshwane (Pretoria)",       city_name := "Tshwane"]
DT[city_name == "eThekwini (Durban)",       city_name := "Durban"]
DT[city_name == "NM Bay (Port Elizabeth)",  city_name := "NM Bay"]

# Colors & legend labels
col_pop <- "#0B3C6F"
col_ntl <- "#F2994A"
labs_legend <- c("Population Density", "Night-Time Light Intensity")

# Sizes
base_sz   <- 12
axis_sz   <- 11
title_sz  <- 13
strip_sz  <- 11
legend_sz <- 11

# Theme
facet_theme <- theme_classic(base_size = base_sz) +
  theme(
    legend.position      = "bottom",
    legend.box           = "vertical",
    legend.box.spacing   = unit(0.3, "lines"),
    legend.text          = element_text(size = legend_sz),
    axis.text            = element_text(size = axis_sz),
    axis.title.x         = element_text(size = title_sz, margin = margin(t = 8)),
    axis.title.y         = element_text(size = title_sz, margin = margin(r = 8)),
    strip.placement      = "outside",
    strip.background     = element_blank(),
    strip.text           = element_text(size = strip_sz, margin = margin(b = 6)),
    panel.spacing.y      = unit(14, "pt"),
    panel.grid.major.y   = element_line(color = "grey85", linetype = "dotted"),
    panel.grid.minor     = element_blank(),
    plot.margin          = margin(t = 12, r = 14, b = 28, l = 12)
  )

# =====================================================================
# FIGURE 4 — Relative to 5 km (bounds 5/10/15/20)
# =====================================================================
df4_pop <- DT[bounds %in% c(5,10,15,20),
              .(city_name, bounds_f, metric = "Population Density", value = rel_pop_d_5)]
df4_ntl <- DT[bounds %in% c(5,10,15,20),
              .(city_name, bounds_f, metric = "Night-Time Light Intensity", value = rel_ntl_d_5)]
df4 <- rbind(df4_pop, df4_ntl)

p4 <- ggplot(df4, aes(x = bounds_f, y = value, fill = metric)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  facet_wrap(~ city_name, ncol = 6, strip.position = "top") +
  geom_hline(yintercept = 1, linetype = "dotted", color = "grey60") +
  scale_fill_manual(values = c("Population Density" = col_pop,
                               "Night-Time Light Intensity" = col_ntl),
                    name = NULL, breaks = c("Population Density", "Night-Time Light Intensity"),
                    labels = c("Population Density", "Night-Time Light Intensity"),
                    guide = guide_legend(direction = "horizontal", nrow = 1,
                                         keywidth = 3.0, keyheight = 0.45,
                                         label.position = "bottom")) +
  labs(x = "Distance from Business District (km.)",
       y = "Density / Density within 5 km of Business District") +
  scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05))) +
  facet_theme

ggsave(file.path(GRAPH_DIR, "Figure_4.png"), p4, width = 7, height = 8, dpi = 300, bg = "white")

# =====================================================================
# FIGURE A1 — Deviation from global average decline (bounds 10/15/20)
# =====================================================================
dfa1_pop <- DT[bounds %in% c(10,15,20),
               .(city_name, bounds_f = factor(bounds, levels = c(10,15,20)),
                 metric = "Population Density", value = diff_global_pop)]
dfa1_ntl <- DT[bounds %in% c(10,15,20),
               .(city_name, bounds_f = factor(bounds, levels = c(10,15,20)),
                 metric = "Night-Time Light Intensity", value = diff_global_ntl)]
dfa1 <- rbind(dfa1_pop, dfa1_ntl)

pA1 <- ggplot(dfa1, aes(x = bounds_f, y = value, fill = metric)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  facet_wrap(~ city_name, ncol = 6, strip.position = "top") +
  geom_hline(yintercept = c(-25, 50), linetype = "dotted", color = "grey60") +
  scale_fill_manual(values = c("Population Density" = col_pop,
                               "Night-Time Light Intensity" = col_ntl),
                    name = NULL, breaks = c("Population Density", "Night-Time Light Intensity"),
                    labels = c("Population Density", "Night-Time Light Intensity"),
                    guide = guide_legend(direction = "horizontal", nrow = 1,
                                         keywidth = 3.0, keyheight = 0.45,
                                         label.position = "bottom")) +
  labs(x = "Distance from Commercial Cluster (km.)",
       y = "Deviation from Global Avg.: Decline from Commercial Cluster (%)") +
  facet_theme

ggsave(file.path(GRAPH_DIR, "Figure_A1.png"), pA1, width = 7, height = 8, dpi = 300, bg = "white")
message("✅ Saved: ", normalizePath(GRAPH_DIR, winslash = "/"))



